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ABSTRACT 

A  generalized  descent  algorithm  theory  is  developed  for  uncon¬ 
strained  minimization  problems.  Here  a  descent  algorithm  is  defined 
as  a  computational  procedure  where  at  each  iteration  a  descent  direc¬ 
tion  is  determined  and  a  single-dimensional  search  is  made  for  the 
minimum  in  the  descen-  direction.  The  theory  is  shown  to  be  a  gener¬ 
alization  of  the  three  most  common  descent  algorithm.’;  gradient,  con¬ 
jugate  gradient, and  Fletcher-Pcwell. 

Execution  of  the  single-dimensional  search  can  be  computa¬ 
tionally  time  consuming.  Two  additional  algorithms  are  presented  which 
reduce  or  eliminate  single -dimensional  search  time.  The  first  is  a 
modification  of  Davidon's  Variance  Algorithm  and  requires  a  minimal 
8 ingle -dimensional  search.  The  second  is  a  direct  method  for  minimiz¬ 
ing  a  special  class  of  quadratic  functions  of  the  form 
1/2 j |xj |2  +  1/2  k(a  -  m'x)2. 
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SYMBOLS 

"*  belongs  Co,  Is  a  member  of;  x£X:  x  is  a  member  of  the 
set  X 

approaches,  converges  to 

+  monotonically  decreasing;  f (xn) +  L :  f(xn)  monotonically 

decreases  to  the  value  L 

[y,zj  value  of  the  continuous  linear  functional  y  operating 
on  the  vector  z.  In  a  Hilbert  space  this  is  the  inner 
product  of  the  vector  y  and. a.  Also,  a  closed  interval 
on  the  real  line,  alternate  to  parenthesis,  and  references 
context  will  make  usage  clear. 

| | x| [  The  norm  of  the  vector  x 

y  The  union  of  sets 

C  is  a  subset  of  (or  equal  to) 

'  matrix  or  vector  transpose  except  in  -Section  XI  where 

f'  is  tne  derivative  of  f. 

Rn  _  iclidean  N-space 

second  differential  operator 

Wij  element  of  the  i£ll  row  and  jit  column  of  tk  >.  matrix  M 
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SECTION  I 
INTRODUCTION 


1.  BACKGROUND 

In  virtually  all  fields  of  the  physical  sciences>and  particularly 
in  engineering;,  the  digital  computer  is  the  principal  tool  used  in  the 
solution  of  complex  problems.  The  speed  and  flexibility  of  the  com¬ 
puter  has  In  many  cases  changed  the  nature  of  the  problems  that  can  be 
solved,  i.e.,  the  solution  must  not  only  meet  specific  constraints  but 
must  also  be  the  best  or  optimal  in  some  specified  sense.  There  are 
three  interrelated  tasks  in  the  formulation  of  such  an  optimization 
problem. 

First,  the  physical  system  or  process  must  be  described  mathemati- 
ec.ly  or  modeled  in  terms  appropriate  for  computation.  Second,  the 
measure  of  goodness,  generally  referred  to  as  the  cost  function,  penal¬ 
ty  function  or  payoff  function  must  be  defined  to  adequately  describe 
hoi,  one  solution  compares  to  another.  Finally,  computational  methods 
must  be  applied  to  find  a  solution  which  satisfies  ti.c  mathematical 
model  and  cost  function  in  such  a  way  so  as  to  extract  the  best  or 
optimum  solution. 

Generally  such  a  problem  can  be  cast  into  a  constrained  optimiza¬ 
tion  oroblem  w,uch  as:  Find  the  solution  x  which  minimizes  the  cost 
function  f(x)  while  satisfying  specific  constraints  described  by 
g(x)=0.  Often  the  problem  can  be  simplified,  conceptually,  by  adjoin¬ 
ing  the  constraints  to  the  cost  function  through  the  use  of  Lagrange 
multipliers.  Thus  the  constrained  optimization  problem  is  converted 
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to  the  following  uncons trained  problem:  Find  the  solution  (x,A)  which 
extremizes  the  cost  function  F(x,A)=f (x)+Ag(x).  The  existence  of  such 
Lagrange  multipliers  is  a  subject  in  itself. 

Another  method  of  solving  the  constrained  optimization  problem  is 
to  restrict  the  problem  to  a  subspace,  an  approximation  to  the  con¬ 
straints  for  example,  and  considering  a  related  unconstrained  problem 
as  an  intermediate  step  in  obtaining  the  solution  to  the  constrained 
problem.  Since  the  simplified  unconstrained  problem  may  have  to  be 
solved  many  times  in  order  to  obtain  the  solution  of  constrained  prob¬ 
lem,  an  efficient  method  of  solving  the  unconstrained  problem  is  essen¬ 
tial.  Finally,  the  solution  of  the  unconstrained  problem  is  often  of 
interest  in  itself. 

The  subject  of  this  thesis  is  the  computational  methods  which  may 
be  used  to  arrive  at  a  minimizing  solution  to  the  unconstrained  minimi¬ 
zation  problem.  It  is  tacitly  assumed  that  any  constraints  are  ac¬ 
counted  for  through  the  use  of  Lagrange  multipliers  or  other  valid 
techniques,  such  as  penalty  factors. 

"  OUTLINE  AND  PREVIEW  OF  SECTIONS 

For  functions  which  have  a  continuous  first  derivative  the  most 
common  methods  used  to  minimize  the  function,  i.e.,  obtain  the  solution 
to  the  unconstrained  minimization  problem,  are  the  gradient,  conjugate 
gradient, and  Fleteher-Powell  algorithms.  These  algorithms  are  reviewed 
briefly  to  illustrate  certain  common  elements.  Here  it  is  assumed  the 
function  to  be  minimized  is  f  which  is  defined  for  each  x  in  some  space 
X.  Further,  assume  the  gradient,  g,  of  f  at  x  also  exists: 
g(x)=grad  f(x).  For  each  algorithm  only  the  initializations  required 
and  the  recursive  equations  are  given.  Convergence  criteria  or  tests 
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for  convergence  although  iaportant  In  computational  applications  are 
omitted  here  in  order  to  emphasize  those  properties  which  these  algorithms 
have  in  coona. 

Gradient  Algorithm: 

Initially:  choose  an  arbitrary  xa 

Iteratively:  set  Sq  =  -g(x^) 

choose  cp'Ojj  to  aininize  fCx^+asn) 
set  x^^  +  ansn . 

Conjugate  Gradient  Algor! tha:  (Reference  5) 

Initially:  choose  an  arbitrary  x,, 

set  s0  “  -gC^) 

Iteratively:  choose  a  =  to  ainiaize  f Cxq+oSjj) 
set  ^n+l  *  *n  *  °nsn 

„  _  |}g(xirf-l)||2 

6»  ilg<vlF 

8n+l  “  ~g(*n+l)  +  Vn- 
Fletcher-Powell  Algorithm:  (Reference  U) 

Initially:  choose  an  arbitrary  xQ 

set  Hn  -  I 

Iteratively:  sn  -  -I^gCxu) 

choose  a  =  ctjj  to  minimize  f(xn+asn) 
set  xn+1  -  xn  + 

°n  -  *n+l  ~  xn 

yn=8K+i>  "  s(xn) 

Hn+l=Hn  “  ^r^n^  +  gnqn 

yA?n  °nyn 

where  the  prime  (')  denotes  transpose. 
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Each  cf  these  algorithms  generates  a  search  direction,  sQ,  for 
which  the  function,  initially  at  least,  tends  to  decrease,  i.e.,  for 
which  g>(xa)sQ<0.  A  single-dinensional  search  is  then  conducted  to 
obtain  the  ainisua  of  f  in  the  direction  Sq  from  the  current  point  Xq. 
The  location  of  the  ainiaua  of  the  single-diaensioaal  search  is  chosen 
as  the  next  iteration  point,  The  differences  between  the  algo- 

rithas  are  in  the  aethod  used  to  generate  the  search  directions  sQ. 
These  algoritfaas  and  others  which  generate  a  descent  direction  and  in¬ 
corporate  a  single -dimensional  search  will  be  collectively  classed  as 
descent  algorithms. 

In  Section  II  the  proof  of  a  theorem  which  is  a  generalization  of 
descent  algorithas  is  presented.  Specific  applications  to  the  gradient 
F'etcher-Powell, and  conjugate  gradient  algoritha  are  given  at  the  end 
of  the  chapter. 

Uezt,  consider  the  problea  of  ainiaizing  the  quadratic  function 
f (x)  *  f  +  a'x  +  l/2x’Gx.  The  gradient  of  f  at  x  is  given  by 
g(x)  *  a  +  Gx.  Let  h  =  -G-1g(x),  assuming  G-1  exists,  then 
g(xtfi)  **  a  +  G(x+h) 

■  a  +  Gx  +  Gh 
*  g(x)  -  g(x) 

-  0 

that  is,  x*  =  x+h  satisfies  g(x*)  =  0,  the  necessary  condition  for  f 
to  have  a  minimum  at  x*.  Note  that  G  is  the  second  derivative  of  f. 
Clearly  for  a  quadratic  function,  knowledge  of  G  or  the  second  deriva¬ 
tive,  or  better  G-1,  greatly  simplifies  the  problem  of  finding  the  mini 
mum  of  the  function.  Since  many  functions  can  be  approximated  by  a 
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quadratic  in  same  neighborhood  of  a  (local)  minimum,  information  about 
the  second  derivative,  or  its  inverse,  should  enhance  the  ability  to 
arrive  at  a  solution  to  the  general  unconstrained  minimization  p  rob  lea. 

The  hewt on-*apfcsoa  algoritfan  is  a  Method  of  function  Minimization 
which  utilizes  the  inverse  of  the  second  derivative.  In  addition  to 
the  computational  difficulties  of  obtaining  the  second  derivative,  this 
aethod  requires  an  initial  estimate  sufficiently  close  to  the  final 
solution  before  convergence  is  guaranteed.  Because  of  these  diffi¬ 
culties,  several  algorithms  termed  quasi— Sevtoo  methods  by  Powell  (Befereict  12 
have  been  constructed  which  iteratively  estimate  the  inverse  of  the 
second  derivative.  The  best  known  of  these  is  the  aethod  of  Fletcber- 
Povell.  Another  more  recent  aethod  of  this  type  is  Davidcn’s  Variance 
Algorithm,  sot  to  be  confused  vitfc  Devi  don’s  Variable  Metric  Algorithm 
which  was  the  predecessor  to  Fletcher-Fovell’s  aethod.  A  new  deriva¬ 
tion  of  Dav.  don’s  Variance  Algorithm  is  presented  in  Section  III. 

Davidon's  algorithn  suffers  some  difficulties  and  Section  III  concludes 
with  a  modified  version  of  the  algorithm  which,  although  somewhat  sore 
complex,  circumvents  one  major  difficulty. 

Whenever  a  method  is  available  for  obtaining  the  inverse  of  the 

second  derivative,  particularly  for  a  quadratic  function,  the  minimi;.- 

ing  solution  can  be  obtained  directly.  In  Section  IV  a  method  for  the 

direct  solution  of  a  special  class  of  quadratic  ninimization  problems 

is  presented.  The  procedure  is  based  on  the  Rank-One  aethod  of  matrix 

inversion.  The  algorithm  contains  a  necessary  and  sufficient  test  for 

the  existence  of  an  extremum  and  a  sufficient  test  that  the  extremum  be 

a  minimum.  The  speci.->l  class  of  problems  to  which  this  method  applies 
are  generalizations  of  the  following  form:  f(x)  =  1/2] | x | j  2  +  k(a-m'x)2. 
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In  this  sect. roc  the  proof  of  a  basic  thecren  on  descent  algorithms 
is  presented  followed  by  applications  to  several  f anil far  descent  algo- 
rit’as. 

1.  SDEfiEIOiS  ASD  5£?'LSniGSS 

Suppose  X  denotes  a  real  no  rased  linear  space  and  f  a  real— valued 
function  defined  an  X.  For  an  arbitrary  point  ^  of  X,  denote  by  S  the 
“level  set"  of  f  at  xq3  i.e. ,  S  *  {x  :  f  (x) <f  (r,,) } .  The  Frechet  deriv¬ 
ative  of  f  at  x  will  be  denoted  f'(x)  and  if  x*£X*,  the  topological 
dual  of  X,  the  value  of  x*  at  x  will  be  written  [x*,  x] . 

let  ♦  denote  a  bounded  nap  fron  S  to  X  satisfying: 

(i)  [f '(x),  ♦(x)j  >  0  for  all  xGS,  and 

(ii)  given  an  e>0  there  exists  6>0  such  that 
[f*(x),  ♦(x)]<5  inplies  ||f*(x)||<e. 

Observe  for  later  reference  that  condition  (ii)  inplies  that 
f'(x)  =  0  whenever  [ f 1 (x) ,  +(x)J  *  0  for  if  there  exists  an  x^  such  that 
If'(xj),  $(xj))  *  0  but  f*(xj)  i  0  set  e  *  1/2 1  }f  *(xj)|  I>0  then 
If'(xj),  ♦(Xj)]  =  0<6  for  all  6  while  J|f,(x1)||~  2e>e  contrary  to  (ii). 
Condition  (ii)  also  inplies  [f'(x),  $(x)]  is  bounded  away  fron  zero 
whenever  f’(x)  is  bounded  a way  froa  zero  in  the  following  sense.  If 
for  {xq^JCS,  |  Jf  * (x,^) |  |>e  for  scae  c> 0;  then  there  exists  a 
subsequence  of  {x^l  and  $>0  suc^  that  [fCx-^),  ^(Xj.^)]^. 

in  the  theoren  which  follows,  -$(x)  serves  to  define  the  descent 
direction. 
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J&ry  of  tise  iiieas  is  tie  following  theorem,  sere  stimnlated 
by  two  risers  by  A-  A-  Goldstein  (References  6,  7).  la 

particular  the  definition  of  $  given  above  acd  tee  form  of  the  coo- 
clusices  of  t be  theor eat  are  identical  to  those  of  Goldstein.  The  hy¬ 
potheses  of  the  theorem  are  changed  to  specialize  to  the  case  of  a 
single-diaeasional  search  at  each  iteration-  Thus  the  proof  of  part 

(a)  of  the  theorem  is  changed.  The  proofs  of  parts  (b)  and  (c)  follow 
Goldstein.  The  following  additional  remark  cm  the  Hypotheses  is  in 
order. 

In  the  current  setting  where  I  is  a  noraaed  linear  space,  the 
assumption  that  f*  is  uniformly  continuous  on  S  may  be  replaced  with 
the  equivalent  conditions  that  f  is  uniformly  differentiable  and  that 
f  *  is  bounded  ca  S  (Reference  li,  p.  ^5)- 

Theoren  I 

Assune  S  is  bounded  in  X,  f  is  bounded  below  on  S  and  the  Frechet- 
derivative  f*  of  f  exists,  is  uniformly  continuous  on  S  and  bounded  on  S . 

Set  x^p  «  Xjj  when  (f '(xq)  ,  $0^)1  *  0,  otherwise  choose  P”Pn  to 
■ininize  {fCx^,  -  p$(*n5)  :p- 0}  and  set  x^  =  x,j  -  pn^(x„).  Then 

(a)  f (Xq)  *L,  f’C^+O, 

(b)  if  {Xjj}  has  cluster  points,  every  cluster  point  z 
satisfies  f(z)=L,  f'(z)=0. 

(c)  If  f*  has  finitely  many  zeros  on  S,  S  is  compact,  and 
llxn41  ~  xnM’*^»  then  the  sequence  (xn)  converges. 

Proof 

(a)  For  xfS  and  f'(x)  i  0,  (f'(x),  $(x)}>0.  In  general,  since  f 


is  differentiable, 


APPEL—' I&-T-— 77 

lf(x),  h]  »  11*  |<f (x+th)  -  f£x)) 
t-*0  c 

»  li*  l(f(*+th)  -  fix)) 
t+0~  t 

•  11*  hf(x-e&)  -  f(x)), 

p*0+  ”** 

so  in  parti.ftu.Sar, 

0  <(f’<x),  Ux) ]  *  14*  %(x-p#(x))  -  f(x)). 

p*0* 

By  the  unifom  differential lity  of  f  there  exists  a  PQ>0,  independent 
of  x,  such  that 

0*-^Hr£(x~p  f(x)>  -  fix)) 

~*s, 

so, 

0>f(x-Po^(x))  -  f(x) 
or 

fCx-pe*(x))<f(x). 

In  particular  for  x  =  x^S,  fCx,.  -  p0^(rn»<f (xQ)<f  (xo>. 

Also  for  f'(x)  /  0,  If’(x),  $(x)J  f  0,  hence  |j$(x)||>0.  Since  S 
is  assumed  bounded,  there  exists  a  Pa>0  such  that  for  all  P>pa, 
x  -  p*(x)fj£s.  For  on  the  contrary  assuoption,  for  every  K,  nc-  natter 
how  large,  there  exists  a  pjj>N  such  that  y  *  x  -  p^$(x)£s.  Then 
1 1  J~x  1 1  *  PN!  U(x)|  I>Nj  j^(x)  j|  is  unbounded  which  contradi  :ts  the 
boundedness  of  S. 

Now  for  Xjjf S  such  that  f’(xn)  i  0,  F(p)  =  fQ 4,  -  P$C*n))  ~  f (*„) 
has  a  minimum  at  Pn£[0,  ca]  since  F  is  continuous  and  {0,  oa]  is  compact. 
Furthermore,  pn  +  0  since  F(p0)  =  f (x,,  -p0<K*n)}  -  fCxn)<F(0}.  Thus 
the  sequences  {pR}  and  {xjj}  are  well  defined  and  f (xn+p  3  fCx„  - 
pn$(xn))<f (xn)  which  implies  the  sequence  {x^CS,  and  the  sequence 
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{f  (Xjj)}  is  strictly  decreasing-  Since  f  is  assumed  bounded  below 
f(*a>  + I~ 

To  shew  f'txjj)  -»0,  we  suppose  due  contrary,  toes  there  exists  a 
sobseamce  {x^}C  J^}  such  that  {'(x^)  is  boaoced  away  fro*  zero, 
lAich  inplies  that  a  subsequence  of  [f'(Xj^),  *(*0^)1  is 
boosded  away  fn»  zero,  Without  less  of  generality  denote  this 
subsequence  by  {1^}.  Then  there  exists  an  eo>0  such  that 

*(*nk)]  *  ^  (f(xnk  -  P+Cx^))  -  ffx^}).  Since  f 

is  unifonely  differentiable,  there  exists  a  p-D>0,  independent  of  x 
(l.e.,  x,,^),  such  that 

T<^b  (Hx°k  ‘  pfc*(*ok))  ‘  f(*ak))» 

>  f {x,^  -  ^(x^))  -  fCx,^), 


or  f(xnk  *  <  f(*nk>  '  f|f°  • 

Bien  f(xuk+i;  “  f(X“k  ~  Pb*(xnk))  <  f(x°k)  "  * 

Since  f  (a^)  *  L, 

fCx^)  <  f(xnk.1)<*“<f(xnk  l+1)  <  f(xD] ,  ^  -  Vo 
t-.at  is. 


£<V  <  £(*»k-i>  . 


<  £:*nk.2>  -  2  . 


£  <W  -  3 


<  f  (s^  >  -  k  V»  . 

**o  2 

This  contradicts  the  assumption  that  f  is  bounded  below,  hence  f'Cx^-K). 
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(fa)  If  *  is  a  cluster  point  of  {xffl}  there  exists  a  subsequence 
{xB^)C{xBJ  such  that  x^^z  (in  non).  Since  f  and  f*  are  continuous, 
f(x„H  and  it  follows  that  f(z)  =  L  and  f'ii)  *  0. 

(c)  Since  f*(z)  *  0  for  every  cluster  point  z  of  {Xg},  the  mmfeer  of 
roots  of  f  *  on  S  is  to  or  greater  than  the  number  of  duster 

points  of  fe  )-  If  f *  has  a  unique  root  z  on  S  then  {x,}  converges  to 
it;  for  otherwise,  since  S  is  assumed  compact,  there  exists  at  least 
one  duster  point  z^  of  {xg}  in  S  where  f'(zj)  0  (by  (b)).  If 
tj  y  z,  the  root  of  f*  is  not  unique.  If  the  umber  of  roots  of  f  ’  on 
S  is  finite  then  the  number  of  cluster  points  of  {Xg}  is  finite  also. 

Let  zj,  i  *  1,2,..., k  be  the  duster  points  of  {xnJ,  let 
e  «  nln{||z£  -  zj||:  i  /  j,  i,j  *  1,2, ...k),  let  SUj.t/S)  denote  the 

open  sphere  of  radius  e/3  centered  at  zA.  Since  S  is  assuned  compact 

k 

the  set  tig}  -y  S(z,,£/3)  contains  a  finite  nuaber  of  poiuts,  say  a. 
i-1 

Since  11*^1  ~  Xj,|  |-*0,  by  ass  imp  t  ion,  there  exists  an  H  such  that 

ll*p*l  "  xpU<  c/3m  for  a11  P5** 

Bow,  since  the  zif  id,... k  are  duster  points  of  Ixq),  there  are 
■eabers  of  {xn>  in  each  S(zife/3)  for  which  u>F  Therefore  fix  n>N  such 
that  xnCS(zi,  e/3),  for  sone  fixed  i,  and  xttrl^  tf(z£,e/3).  Let  xq 
be  thi  next  member  of  the  sequence  }xn|  (i.e.,  q>n)  such  that  xq£_S (zy 
e/3),  for  some  j  /  i  and  xq_^^S(zj  ,e/3). 

Since  xnCS(.zi,c/3)  and  xqCS(zj,  e,  3) ,  i  /  j ,  |  |xn  -  r  J|>e/3. 

On  the  other  hand,  since  n,  q>N,  we  have 


l*q  -  xn I  1=1  !  xq  -  xq-l!l  +  Hxq-1  "  *q-2 • 1  +  ’** 

+  Nxn+2  ‘  *n+lH  +  I|xn+\  ~  *JI 

■“i1 

p=n  p=n  J  J 
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Then  ^  <j  jx^  -  xj  |  <  (q-n)  ^  cr  e<f-a. 

low  the  set  of  points  oeloot  r°  the  set 

k 

{*_}  -  u  S(z, ,  t/3).  lost  there  are  (q-1)  -  (nfi)  +  I  *  q— a-1  points 
1-1 

tc  the  set  tifl| ***,  x^_j),  aa£  by  supposition  n  points  in  the  set 
k 

W>  -  u  S(z- ,  c/3).  Thus  On  the  other  half  we  have  just 

i-1 

shown  q— n>n,  thus  s£Lq— m-l>sr-l,  that  is  to  say  q-o-1  =  a.  Heuristically, 

k 

all  the  points  of  the  set  (x_)  -  y  S(z»,  c/3)  have  been  accounted  for, 

i-1 

or  “used  up".  But  by  the  same  argtssent,  there  exists  an  n*>q  and  a 
q*>n'  such  that  xQ» CS(zj ,  e/3)  and  £s(zj ,  e/3),  for  the  j  de¬ 

fined  by  q  above  (i.e.,  for  which  x^  £S(Zj,  e/3 )),  and  x^.Csfz^,  e/3), 
1  j,  and  ^ S (z^ ,  e/3).  Then  there  are  again  q’-n*-l  *  a  points 
in  the  set  (x^*.^,  xn'+2*  *  *  * »  xq*-l^  also  belong  to  the  set 

k 

{x  }  -  u  S(z, ,  e/3).  There  now  are  2a  points  in  this  set  which  by 
i-1 

supposition  contained  only  a.  This  contradiction  persists  unless  we 
suppose  {Xq}  has  a  unique  cluster  point, 

Coroaents  on  the  Basic  Theorea 

The  first  conclusion  of  the  theorem,  f(xn)+L,  can  be  obtained 
with  the  weaker  assumption  that  f  is  differentiable  and  without  the 
assunption  that  $  is  bounded,  in  which  case  both  pQ  and  pa  are  dependent 
on  x  (or  xq).  However,  both  uniform  differentiability  of  f,  which  fol¬ 
lows  from  the  uniform  continuity  of  f,  and  boundedness  of  $  were  used 
to  assure  the  existence  of  a  which  uniformly  bounds  the  differential 
away  from  zero  in  the  proof  leading  to  f'(xn)-^0. 


A?75Lr-33-T2-T7 

If  X  is  finite  dlansioul  tbcc  S  being  closed  and  bowled  Is 
necessarily  rnrparr  sad  Ixfl}^S  hit s  cluster  points.  Bence  conclusion 
(b)  of  the  theorem  applies. 

3-  d??IIC2SIC3  3?  TH2  -ASIC  TES3H2X  SO  CMOS  S2SC23P3F  fLOOSTSEJS 

In  the  applications  whi ch  follow,  the  boundedness  of  $  will  be  de¬ 
rived  f  ran  the  boundUess  of  the  Frechet  derivative  f  *  cf  f  on  S.  In 
particular,  assone  now  that  X  1s  a  Hilbert  space  so  that  f '(r)  be 
represented  by  its  gradient  Vf  00  in  X.  Then  also  since  in  a  Hilbert 
space  every  bounded  set  is  weakly  compact,  the  boundedness  of  f  *  on  S 
follows  iron  its  continuity  on  S(Eeftreice  .3,  p.  19). 

Corollary  1 

Let  Q  be  a  positive  definite  continuous  linear  operator  on  X, 
let  #00  «  Q7f(x),  then  Theoren  I  applies. 

Proof: 

Since  f  has  the  required  properties,  all  that  remains  is  to 
show  +  is  bounded  and  satisfies  conditions  (i)  and  (ii) . 

By  assumption  f  *  is  uniformly  continuous  and  bounded  on  S.  This 
together  with  Q  being  a  continuous  linear  operator  on  X  implies  4  *  QVf 
is  bounded  on  S. 

Q  positive  definite  implies  there  exists  a  n>0  such  that 

■IMI2  <U»  QD  *---  a11  2  in  x* 

hence 

If*  00,  $00)  =  If' 00,  Qf ’  (x)]>,  n|  ’  (x)|  j2  >p  and  condition  CD 
is  satisfied. 

To  show  condition  (ii)  assume  an  e>0  and  choose  c'=me2.  Then 
(f'(x),  ♦(x)J<  6=me2  implies  m| jvf (x) | J 2^{Vf (x) ,  QVf(x)J<  6=me2  or 
j ;Vf (x) | |<e. 


■pv 
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In  order  Co  demonstrate  that  Che  theory  developed  tints  far  is  a 
general izat ion  of  descent  algorithms,  the  Iheoret  and  Corollary  will  be 
applied  to  two  comipa  algorithms ;  the  gradient  and  the  Fletcher-Powell 
methods. 

Application  to  Gradient  Algorithms 

Let  Q  in  Corollary  1  be  the  identity  operator,  i.e. , 

♦GO  2  7f(x),  this  is  the  usual  gradient  algorithm. 

lae  operator  Q  need  not  reaain  fixed  so  long  as  it  is  uniformly 
positive  definite  on  S.  That  is  to  say  there  exists  a  constant  *>0, 
independent  of  x,  such  that  If*G=)»  #GOJ  *  If'(x),  QOOVf  GOJ>p|  |vf  (x)  j  |2 
for  all  x  in  S- 

Aoplic^ioc  to  Fletcher-Powell 

Let  X  be  finite  disensi<aal,let  Q  *  Q(k)  *  as  defined 
by  R.  Fletcher  and  M.J.D.  Powell  in  Reference  h.  Using  the 
algorithn  of  Fletcher-Powell  it  is  easy  to  prove  remains  positive 
definite  at  each  iteration  unless  the  algorithm  terminates  in  a  finite 
number  of  steps.  Thea  ♦(Xq)  =  H^f (Xg)  satisfies  condition  (i): 

{f'CXjj),  ♦(xn)J>0.  Sow  if  the  are  uniformly  positive  definite,  or 
if  there  exist  constants  m>0  and  p>0,  independent  of  k,  such  that 
Jx,  Hjjx]^.  nj  |xj  p,  then  condition  (ii)  is  satisfied.  For  given  an 
e>0,  choose  6=a  ep,  then  m| |f ' (x) [ |p^[f '  (x) ,  H'Cf'(x)J<  6=meP  implies 
| j f  * (x) 1 1 <e .  Then  the  theorem  applies  to  the  Fletcher-Powell  algorithm. 
Furthermore,  since  S  is  closed  and  bounded,  it  is  compact,  and  the 
sequence  {x^}  necessarily  has  cluster  points  so  (b)  of  the  theorem  ap¬ 
plies  ithout  further  assumption. 

In  the  above  application  to  Fleccher-Powell 's  method  and  in  the 
following  corollary,  the  function  $  is  not  uniquely  defined  as  a 
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function  fro*  S  into  X.  lather,  )  is  defined  only  on  the  s-^quence  {^J 
and  its  value  at  an  arbitrary  x  in  S  depends  upon  which  sequence  x  be¬ 
longs  to.  However,  the  properties  of  #  used  in  the  proof  of  Theorem  I 
did  not  depend  on  $  being  defined  anywhere  except  at  points  of  the  se¬ 
quence  {Xjj}.  That  is  to  say  the  conditions  (i)  and  (ii)  need  hold  only 
at  the  points  of  the  sequence  {xq}.  Within  this  context  the  search 
direction  will  continue  to  be  denoted  by 
Corollary  II 
Let  ♦(xQ)  -  and 

♦(*a)  *  VfCxjj)  +  for  n  -  1,2... 

where  |Kn|<y<l,  then  Theoren  I  applies. 

Proof: 

Clearly  conditions  (i)  and  (ii)  are  satisfied  at  xQ  for 
C f ' (x0) ,  +(x0)  )  =  |jvf(x0)jj2  >p,  and  given  e>0  choose  6  «  e2. 

Claim  ♦(*,!)]  *  0  for  n  *  0,1,....  Since  pQ  is 

chosen  to  minimize  (F(p)  =  f  (xj,  -  pfCx^)  -  f^):  p>0),  F' (pn)  *  0. 

But  F'(p)  »  -  [£*  (Xp  -  p+Cx^)),  ♦(Xjj))  and  at  p»pQ,  x,j  -  p^Cx^  -  x^. 
Hence  F*(pQ)  =  -  (f'Cx^j),  ^(Xjj)]  -  0. 

Then  for  n  =  1,2,. . . 

(f'(xn),  3  *  [f^Xjj),  VfCxjj)]  +  Kjjtf'Cxjj),  (Kxjj.j)] 

-  (Vf(xn),  Vf^)] 

°  ||Vf(Xn)||2>0. 

2 

Hence  again,  given  e>0  choose  6  -  c  then  [f'Cxjj),  ‘Kxjj)]  <  6  implies 
|jvf(xn)||  =  ||f,(xn)||  <  c.  Thus  conditions  (i)  and  (ii)  are  satisfied. 

To  show  boundedness  of  the  «f>  (xjj)  consider 
lj^(xn)||2  =  (^(Xjj),  $(*n)]  =  [Vf(xn),  $(*„>]  +  KnU^-l*’ 

“  llvf(Xr.)||2  +  Knl,Kxn-l)»  4’(xn)],  and 

Ih 
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[♦(*n).  ♦Gfe-1)]  =  [Vf(xn>,  <t'(xn_1)]  +  K^Hx^),  ♦(xn_1)J 

IN*n-1>ll2- 

Thus,  |U(xn)||2  =  ||7f(xn)||2  +  K2||Kxn_1)||2. 

Applying  this  relation  recursively  yields 

IN(xn)||2  =  |  IvfCxn)!  |2  +  K2(||Vf(xn_1)||2  +  K2_1||Kxn_2)||2} 

“  IIVfCxH)!!2  +  K2||7fCxn.1)||2  +  K^K2_1||Hxn_2)||2 

•  •  • 

I I+c^s,)  1  !2  *  llmx^ll2  +  k2! 1 2f o^i) 1 1 2  +  ^g_x|  | Vf (x^) 1 1 2 

+  •••  +  K^K2^!  •••  K| 1 1 Vf  0c1) 1 1 2 

+  KnKn-l  —  KlKfl  |vf  0co>|  |2- 

Since  f'(x)  is  bounded  for  xCS  let  H  *  sup  ||vf(x)|| 

xCS 

then 

IlKXn)!!2  <  U+  K2  +  K2^  +  •••  +  K2^-!  •••  K^M2. 

Also  K2  ^  r  <  1,  so  the  series 

{!  +  K2  +  K2^  v  —  +  KjK2.!  K2*2} 
converges.  Hence  ||(Kxn)||2  is  bounded  and  therefore  is  bounded. 

For  a  trivial  application  of  Corollary  II,  Kj,  may  be  set  to  zero 
for  all  n,  this  then  generates  the  usual  gradient  algorithm.  The  fol¬ 
lowing  application  is  of  much  more  interest. 

Application  to  Conjugate  Gradient 

Let  Kn  =  IVf  (*n)  ,Vf  (^>3 _  in  Coronary  n. 

!Vf(xn-iX^f (xn-i)J 

This  is  the  Bn-i  of  Fletcher-Reeves  (Reference  5)  for  X  =  Rn,  and  Lasdon, 
Witter  and  Warren  (Reference  10),  for  X  a  function  space. 
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The  condition  IlCjJ  <,  r  <  1  nay  appear  overly  restrictive  particu¬ 
larly  in  the  light  of  the  paper  by  Lasdon,  Hitter  and  Warren  which  also 
is  set  in  a  Hilbert  space  K.  Unfortunately,  the  proof  of  Lasdon,  Mitter 
and  Warren  contains  one  minor  but  signficant  discrepancy.  The  proof  of 
their  theorem  is  reproduced  and  corrected,  to  illustrate  the  simi¬ 
larities  of  the  constraints  which  must  be  imposed  on  or 
Theorem  3  of  Lasdon ,  Mitter  and  Warren 
If:  1.  J(u)  is  bounded  below, 

2.  J(u)  and  g(u)  ■  grad  JCu)  are  continuous, 

3.  D2J(u,h,h)  exists  and  |D2J(u,h,h)|  £  ra||h|pfor  m  >  0  and 

all  u,h  in  H, 

4.  {uk}  has  a  cluster  point  u*, 

then  the  sequence  {u^}  formed  with  arbitrary  uQ  by  applying  the  con¬ 
jugate  gradient  algorithm  to  J(u)  has  the  following  properties: 

1.  lim  J(uk)  *  J(u*), 
k->« 

2.  lim  g(uk)  ■  g(u*)  -  0. 
k-*» 

Before  presenting  the  proof  the  following  three  remarks  are 
pertinent: 

1.  Here  D2J(u,h,k)  is  the  second  differential  of  J  at  u. 

2.  The  form  of  the  algorithm  is  exactly  as  given  in  Section  1.2 

except  here  the  independent  variable  is  u  instead  of  x. 

3.  From  the  above  proof  of  Corollary  II  it  follows  that 

(8k>skJ  “  -!|gkll2  and  fsk,sk+l^  =  0  for  a11  k- 

Corrected  Proof  of  Theorem  3  of  Lasdon,  Mitter  and  Warren 
From  Taylor's  Theorem, 

J(uk  +  ask)  =  J(uk)  +  a[gk,sk]  +  l/2a2D2J(£k(a) ,sk,sk) 
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where  £00  belongs  to  the  line  segsent  joining  uk  and  Uj+ask.  Then, 
using  Isk,gk]  =  -IJgjJl2  and  assumption  3, 

J(uj.+osk)  <  J(ufc)  -o|lgkl|2  +  1/2AJ  |skS|2- 


Since  a?=ak  nininizes  J0»k  +  ask>, 

JO^  +  o^s^)  =  J(uk_j.^)  i  ^ 

<J(uk)-i||gkll2+^li»kll2- 

At  this  point  Lasdon,  Hitter  and  Warren  assume  j |sk[ |  =  ||gk|| 
in  their  proof  of  this  theorem.  This  is  clearly  inconsistent  with 
[sk,gk]  =  — | *gk| I 2  f°r  by  the  Cauchy-Schvarz  inequality 
Usk’8k3l  =  !|skM  !  Uk! !  wb^re  equality  holds  if  and  only  if  sk  is  a 

multiple  of  gk-  Thus  !  1  Sk 1 1 2  =  I  Isk»f>kJ  i  £  1  !skl  I  Lgkllor 
| |gkJ  J  <  ||sk|j.  Specifically,  equality  holds  only  when  k=0 
(i/e. ,  sQ  **  -gD)  or  =  0  in  which  case  the  solution  has  been  achieved. 
Proceeding  more  carefully,  observe  that  sk  =  -gk  +  Pjt-^k-l’ 

!|skll2  *  Isk’sk3  =  ~Isk’gk3  +  ek-lIsk»sk-l3 
“  Ilgji2  +  gk-l^sk’sk-l3 
and  Isk-!*skJ  =  -Isk-l,8k3  +  ek-llsk-l’Rk-l3 


■  ek-l*sk-l’sk-l3’ 

or  |  |sk|  |2  =  |  |gkl  I2  +  ^-i  |  Isk-ilK 
Then 

J(uk+1)  <  j(uk)  -  i| lekl I2  +  ^<l|gkli2  +  6k-l  llsk-lll2> 

-2S|l8kl|2+S°k-i  l|sk-i|i2' 


cewise, 

J (uk)  ~  3(uk_i^ 


2^  ^8k-ll 


_IP2 

2m  k 


-1 


1 1  sk—  2  I 
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hence 


«Vl>  -  JC“k-l)  ‘  Si  I2  +  S6^!  K-2I  l! 


-^M**£Wiviii2- 


1  qZ 


or  in  general. 


V  i  J<V  - 1  s'M2  +  X  s'lM'iM2- 


or 


k-X  0?  k 


J<V  <  J(uQ)  +  I  gllsjl2-  I  ^NgJI2 


i*0 


i-0 


'1-slkll2. 


1-1 


where  zt  -  _I  1 1,  !!'  -fLi  1  Ni-lli2- 

2a  2a 

k 

Since  J(u)  is  bounded  below  lin  £  z^  exists  and  is  finite,  that  is  to 


k-M.  i«i 


say  the  series  Jz^  converges. 

Since  \z±  !g±!  I2  *  ^Ulsi-il!^]  it  follows  that 

(Reference  13)  if  the  series  Jz^  is  absolutely  convergent, 

then  each  series  1 1  1 1 2  and  [s._,  [  | 2  is  convergent,  which 

2m  1  2m 

implies  JlgJI-0  and  &2-il  lsi-il  1 2"*0»  It  also  follows  that 
llsi!l2  “  M g± 1 1 2  +  ihi.ill2^)  from  which  it  follows  that 

Msjl  is  bounded.  On  the  other  hand  if  \z^  is  conditionally  con¬ 
vergent  (convergent  but  not  absolutely  convergent)  then  each  series 

yikll2  and  I5Ll.ll*i  _jjj2  is  divergent.  Of  course  it  still  may 
2m  1  2m ' 

happen  that  ||gjJj-*-C  and  8|„jJ  lsi-i I  I  2"*0  (and  hence  S£  is  bounded). 

One  way  to  assure  £z^  is  absolutely  convergent  is  to  assure 
z^O,  that  is  2ra  z±  =  MgJI2  -  8^1  1  SJ_-1  i  I^P-  But 

Msi 1 1 2  =  |  |e± 1 1 2  +  32_i  lisi-ill2>  hence  2m  zi  =  2 1  Ui 1 1 2  -  l|si||2>o 
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re«|tilres  ||sjl|2  <  ZW^U2  (recall  ISgJ!2  <  IfsJI2).  *cr 
llsjl2-  IlgJI2, 

IKII2-  lUill^eJilsJI2 

*  !lgj|2  (1  +  3  M) 

u  °  ,,*oH2 

“  IIrxI!2  a  + b0) 

since  8o  =  | |g3 | |2/| JgJ |2- 

Similarly, 


1I«2IP-  1U2IP  a  +  vJf^] 

-  I|g2ll2  a  +  aLa  *60)> 

-  Ilg2ll2  (1  +  8!  +  8j80) 


or  in  general 


I  l*kl  I2  *  I  lgfc.1 12  (1  +  Bk  +  Mk-l  +  —  +  h.h-1  ***  8  > : 


Now  if  eA  <  r  <  1/2  for  all  i. 


1  ♦  “i  +  Vk-1 +  -  +  skVi  6oiw<'Firt  *  2 

Then  i  lskl  1 2  <  2 1  |gk|  1 2»  and  zk>0,  and  absolute  convergence  of  £*k 

follows  from  its  convergence.  The  condition  8^  <,  r  <  1/2  is  stronger 

than  the  condition  ] |  <,  r  <  1  assumed  in  Corollary  II. 

N 

Application  to  Conjugate  Gradient  in  R 

The  classic  proof  of  the  conjugate  gradient  algorithm  for  function 
minimization  in  finite  dimensional  spaces  (R1^)  is  based  on  the  Gram- 
Schmidt  Orthogonalization  procedure.  Originally  the  method  was 
developed  for  the  solution  of  systems  of  linear  equations.  The  ex¬ 
tension  to  the  problem  of  minimizing  a  quadratic  function  on  is 
well-known.  In  theory,  the  conjugate  gradient  method  finds  the  minimum 
of  a  quadratic  function  in  at  nioct  N  steps.  However  in  practice,  be- 


iSV5L-'VSr-T2~'T7 


cease  of  roosd-off  error,  X+l  steps  are  etei  to  obtain  rise  "exact"' 
solution. 

Application  of  rise  ooajtrgate  gradient  algorithm  to  rise  problem 
of  minimizing  an  arbitrator  function  in  S?  usually  follows  toe  same 
procedure,  viz. ,  conjugate  directions  are  used  for  Srt  steps  and  if 
the  minimum  has  not  been  obtained,  rise  algorithm  is  re-started  with 
the  last  best  estimate  (^).  In  teems  of  Corollary  II,  this  means 
Zgll  *  *SfJH-2  “  1  ®*  This  condition  may  be  tssed  to  impose  a  some¬ 

what  weaker  but  less  instructive  condition  on  the  ^is  in  Corollary  II, 
specifically 


a  -Sj  +  «£i.1+~  +  *i«i.1~*i>4A 

«  +  «an  +  +  —  +  K2KHrl*'"1^2>  s  *  «<=-  f“ 

each  sub-cycle.  Then  j  ]  ACxfc)  |  j 2  <A  M2  for  all  k  where  J  |vf  Or,  )  |  J  2<.  M2. 

As  a  final  remark  it  should  be  noted  that  in  Corollary  II  the  con¬ 
dition  |kJ  <.r<l  need  not  be  satisfied  for  all  n  but  only  for  all  n 
beyond  some  point-  In  the  co-  jugate  gradient  algorithm 

K_  *  i fo^fon) 1 1  <  r  <  1  implies  the  gradient  not  only  converges  to 

llmxn.1)||2 

zero  but  | } Vf (x„) | J  forms  (eventually)  a  strictly  monotone  decreasing 
sequence. 
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HCTIt  TTITT 

&  aaas-gas  m shjs  a?  HBGrsni  xmxFSEn-gg 

A  particular  property  possessed  by  kotk  tie  aetfeol  of  cnjapte 
fradiats  xd  the  Method  of  Fletdser-fovell  is  that  either  Method  ob¬ 
tains  the  ■<»<—  of  a  positive  definite  quadratic  fom;  viz.,  mini— 

Mizes  f 00  *  f  +  a'z  +  1/2  x'Gx,  xO®,  is*  a  finite  wither  of  steps 
excepting  rocod-off  errors.  This  is  accomplished  in  the  cnajogate 
gradient  algorithm  by  means  of  the  Gram-Schmidt  orthogonal  1  ration 
process  (asffiresce  1).  In  the  method  of  FLetcfcer-Pouell ,  this  Is 

aerwpl ished  by  generating  the  inverse  of  the  Matrix  G,  specifically, 

Bg  ■  C T1.  In  thin  section  another  method  for  generating  G~*  Is  pre¬ 
sented  which  does  not  require  a  single-dimensional  search  for  a  Minima 
as  do  the  methods  of  Fletcber-Powell  and  conjugate  gradients. 

Although  this  author  arrived  at  the  method  independently,  the 
algorithm  is  essentially  the  came  as  Davldon's  variance  algorithm  (Reference  3). 
Although  David  on's  proofs  are  valid,  they  provide  little  insight  for 
the  user  on  how  the  algorithm  is  structured  or  why  the  method  works. 

The  derivation  presented  here  clarifies  the  structure  of  the  algo¬ 
rithm  by  emphasizing  how  the  structure  leads  to  the  desired  properties 
at  each  iteration.  Finally,  a  minor  change  is  incorporated  into 
Davldon's  algorithm  which  circumvents  one  difficulty  which  may  be  en¬ 
countered  when  the  algorithm  is  applied  to  a  computational  problem. 

1.  BASIC  RANK-ONE  METHOD 

In  Reference  15,  Herbert  S.  Wilf  presents  a  netnoa  for  matrix  inversion 

based  on  the  equation  (A  +  uv')  *  =  A  ^  — 

1  +  v’A-1u 

uhere  the  prime  (')  denotes  matrix  (vector)  transpose  and  it  is  assumed 

A"^  exists.  This  equation  is  easily  verified  by  computing 
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(A  +  w*)C A"1  -  *  — >A  I  1  «  I. 

1  ♦  t'l'V 

Sccc  t!u:  A  +  sbt*  differs  fraa  A  by  a  Matrix,  «v*  ,  ef  rad  sue.  litre* 

fore,  tie  abere  issreisioe  terfaai^me  aad  tie  algorithm  developed  below 

sight  veil  be  called  aagfcHjjtae  Methods. 

Folltvijj  tfilf ,  ssspp ose  C  is  written  ia  tie  form 

«  a 

C  *  C0  +  J  «!*!  aad  let  be  tie  partial  s«s  C,  1  C.  +  |  nf»!  aad 
i-1  1  1-1 

Da  denote  the  partial  sob  inverses 

^  ^ 1  -  cco  +  *  o^r1. 

i-l 

n“(Wl'^  +  Vrfti 

«♦*  _1 

Dr*l  -  <Co  + 

“  CCn  +  WW1 


-  D  -  Dnan*-lvgfrlpn  . 

°  1  +  vn+lDnun+l 

-1  M 

It  follows  that  Dl.  =  tC  -  (C  +  j>  u,v. )  -  CT*  providing  all  ccopuca- 

i-1 

tions  *ay  be  carried  out.  Obviously  the  procedure  fails  if  the  denom¬ 
inator  tern  1  +  is  zero.  The  inplication  of  such  a  condition 

can  be  deduced  from  the  results  of  the  following  Lera a. 

Lean  a  I 

Assure  A  is  nonsingular,  then  the  natrix  B  =  A  +  uv'  is  singular 
if  and  only  if  1  +  v'A~^u  =  0. 

Proof 

Necessity :  Assume  B  is  singular,  then  there  exists  an  x  /  0  such 
that  Bx  =  Ax  +  uv'x  =  0.  Then  A~*Bx  =  x  +  A~*u'-,x  =  0  and 
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*  r*x  {1  +  w*A”*w}  *  0.  Sew  r'x  f  0,  for  tie®  8x  *  Ax  *  0 
Wt  A  Is  acaslncnlxr-  Tkerefere,  (1  *  t*&s5  *  0. 

$*£ ficfescy:  Assume  1  +  v*A®  *  SB  ami  set  x  *  A  /  ©  for  cdser- 
wise  o*0  v*A— -  0  which  cfflctradiets  I  ■#■  »*A"^c  *  0.  The® 
lx  *  Ax  +  mr’x 

*  n  +  mr*A~^o 

*  (J  +  *'i  ^a)  «* 

-  0. 

2.  S52T.2CT22  0?  5=2  £02-352  X3DGZ£21fflS  AL/30EI5EK 

An  algorithm  will  be  anstrnctd  for  minimizing  tfce  quadratic 
fo,x  f(x)  *  f 0  +  a*x  +  1/2  x'GjC,  xC^*  and  C  a  positive  definite  sjb- 
aetric  matrix.  The  method  is  based  co  the  above  technique  for  generat¬ 
ing  Cf1.  Let  Lq  *  CT*,  where  Cc  is  an  arbitrary  positive  definite  sym- 
*etric  oatrix  such  as  the  identity  I.  Pick  an  arbitrary  xQ  and  let 
gQ  *  a  +  GxQ  =  grad  f(x0).  As  suae  gc  /  0  for  otherwise  xQ  is  the  required 
solution.  Set  aQ  *  -D^,  Xj  *  xD  +  oQ,  gj  =  grad  fCxj)  =  a  +  Gxj  and 

y0  -  si  -  g0-  Then  si  =  a  +  G<*o  +  °o>  =  «o  +  too  2nd  yc  =  ^o' 
improved  estinate  to  G  is  sought  such  that  Cj0o  *  CcQ  =  yQ,  where 
Cj  has  the  forra  *  C0  +  u^vj.  Since  G  is  symmetric,  being  the 
Hessian  of  f,  the  added  constraint  u^  *  is  imposed  tc  assure  sym¬ 
metry  of  C^-  Then  the  condition  to  be  satisfied  is  CjOD  =  co0o  + 

ului°o  =  yo»  or  U1  =  ~  'ip—  *  But  ’  Co°o  “  ~Co(-Do8c,)  =  «o»  50 

ul°o 

Ui  *  _5i—  =  kgi.  The  value  of  the  scalar  k  is  easily  determined  from 
ul°o 

8i  ?  1 

kg,  =  — ■*—  to  be  k^  =  r~  and  Ci  takes  the  very  simple  form 

1  ui°0  81  0 

g  g1 

Cj  =  CQ  +  lbl  .  Applying  the  rank-one  inversion  formula  to  yields 

g’o 

1  o 
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HmCl 


er 


Vi  c-i 

^ .  <r"  -  1  *  , 

c  , 

I*3£l 

*1°© 

si  Vi  **iDofe.  +  *o> 

-  *iVo  '  «i®0* 


*°  ®1  *  ®0  ”  C  j;  1  2  . 
*iVo 


Bost 


0  g  *’9 

-  olio 


#  ‘iCo  *  *Po*l 


Ju*t  as  C1<a0  -  yQ,  D^y  -  cQ  for  l^y,,  -  D j0  -  DoZlSiDc»o 

*iVo 


“  D0fel  "  V  -  DcSl  -  -»oSo  “  % 

or  simply  by  applying  ^  to  Cj0o  •  yQ. 

Since  Dj  is  the  current  best  estimate  of  C  \  it  is  reasonable  to 
continue  the  process  by  computing  o^  *  ~D1&1»  x2  *  X1  +  °1» 
g2  *  grad  f  (xj)  and  =  §2  ”  ^1  ^  considering  an  improved  estimate 
<2  to  G  obtained  from  C^,  g£  and  y^  in  precisely  the  same  manner.  More 
generally  suppose  xn,Dn,gQ  and  CQ  have  been  obtained-  Set 

°n  * 


*n+l  “  *n  +  °n» 

Sn+!  =  8rad  f<*n+l> 


and 


~  ®n+l  ~  6n 


Go„. 


Define  =  DQ  -  Dngn+lgn+lDn 

^n+l^n^n 


and  Cn+jL  =  Cn  +  gn+lgn+l 


gn+l°n 


Clearly  Cn+1on  =  Cnon  +  gn+1  =  -gn  +  gn+1  =  yn  =  Gon  and 
Dn+iyn  =  Dnyn  "  Dngn+1  =  DnSn+l  "  DnSn  '  Dngn+1  =  °n  =  G_lyn- 
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Ac  iagoircaas.  property  of  tSae  algorithm  is  if  fl^y  *  a  m  C  *y  for 
as  arbitrary  y  tbes  ^afrlT  *  ^s?*  ®>f»  property  is  demonstrated  by 
shovI=s  *  0  as  follow.  Since  g^  -  *»  *  *0^*1  “  **)» 

cs  *  vr  xa  1  c_1teo-i  -  *»)  *  -aa». 

«  "  C_1Csffi4-i  -  «a>  *  0. 

Ihe®  ~  c"1tSafrl  -  8b5  +  ®sA*+l 

*  0>n  '  ^Htn+l  '  «o>. 

«■*  4+lV  *  <*o+l  "  *n),(In  -  C_1>y  *  °* 

That  is  to  say,  if  I>n  agrees  with  C-^  for  sok  vector  y,  risen  so  does 
D«*l* 

Clearly  cannot  be  computed  if  the  ten  gV_jDaya  vanishes. 

If  this  occurs  as  above  because  DQ  agrees  with  C"1  on  yQ,  then 

*  c~\  *  ^(Sn+l  ~  *n> 

«  (T1(a  +  -  a  -  GXjj) 

“  *n*l  *  *n  “  °n  *  "Dn*n*  But  also*  Vn  *  Dn*n+1  "  Vn’  be¬ 
fore  *  0  and  if  Dn  is  non-singular  g^fi  “  0  which  in  turn  implies 

is  the  desired  solution. 

By  Leona  I,  Dul  ^  exists  so  long  as  the  denominator  (1  +  v’A  *u  in 
the  Lexsaa)  does  not  vanish.  In  the  case  of  the  D£s,  ;  exists  so  long 

35  *n+l°n  +  CnDn8n+l  =  S^lVn  *  °-  Furthermore,  can  be  ob¬ 

tained  from  Dn+1  by  applying  the  same  inversion  technique,  hence  the 

denominator  in  the  expression  Cn+^  =  CQ  +  *rH-l*n+l  nust  not  vanish. 

*n+l°n 

There  are  then  two  requirements  for  the  existence  of  Dn+1  =  <£+1  and 

Cn+l  =  D;+l»  naMely=  CM  *  0  and  *n+l°n  *  °- 
In  the  following  it  is  assumed,  for  the  moment,  that  these  require¬ 
ments  are  satisfied. 


2> 


TO-T?— ' V* 


The  tltoritta  ur  fee  (faciallztd  to  one  for  ■tafeiztag  a* 
arbitrary  fsisctioo  f  defined  oa  Jp  as  follows: 

Initially:  qook  an  arbitrary  xQ,  set  D0  m  I 

sad  compare  gQ  *  grad  f  (r^. 

Iteratively:  Sec  ®n  *  -D^  aad  +  ®n- 

Conpate  gn5.1  *  grad  f  (x^+j) » 

***  7»  *  *n*l  *  *n  »* 

D  g  g*  D 

compute  D^j  *  DQ  -  n  tri-1  afl  n  . 

^flVn 

is  an  algorithm  for  alninlzing  an  arbitrary  function,  precautions 
xsst  be  taker-  to  avoid  instances  where  the  deossiattor  sill  vanish. 

In  addition,  for  an  arbitrary  function,  Dq  nay  not  be  positive  definite 
and  hence  aQ  *  _Dn^n  aay  not  a  descent  direction.  Therefore,  addi¬ 
tional  precautions  must  be  incorporated  into  the  algorithm  to  assure 
applicability  to  arbitrary  functions  at  the  expense  of  increased  com¬ 
plexity. 

The  recursive  relation  for  nay  be  rewritten  as  follows: 

Vi  -  ®.  -  "AhCA 

E^“»y"  CAVi 


or 


n+1 


+  (A  -1)  Dn8n4-lg^-lDn 

»  *1  T\  ^ 


gn+lDn8n+l 


where  =  1  -  gn+lDngn+l  . 

ClVn 

This  form  for  Dn+1  is  essentially  the  recursive  relation  of  Davidon's 
Variance  Algorithm  except  tha".  X  is  chosen  to  assure  both  Vi  and 
D~^  (i.e.,  Cn+^)  remain  positive  definite  at  each  iteration,  where  tiie 
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above  relation  for  2^  is  used  whenever  possible.  Specifically,  1Q  is 
chosen  seen  that  for  all  x,  where  0<a<i<g- 

As  Dari  doc’s  Variance  Algorithm  is  presented  in  Heference  3  , 

Vi  (or  ^  in  Daw  idea’s  notation)  is  constructed  from  a  test  value  x* 
for  If  f(x*)  >  fCxjj),  then  ^  is  taken  for  that  is  the 

estimate  x*  is  discarded.  As  pointed  out  by  Davidon  in  a  footnote 
(Seferenee  3,  p-  +06}  the  algorithm  can  become  trapped  in  a  loop.  Indeed, 
although  the  poor  estimate  x*  is  not  used,  the  gradient  of  f  at  x*  is 
used  to  modify  the  estimate,  D^,  of  G T*.  Since  for  most  problems  f(x) 
can  be  computed  nuch  more  rapidly  than  g(x)  *  grad  f(x),  several  test 
computations  of  f(x),  without  the  corresponding  gradient,  can  be  made 
without  undue  increase  in  computation  time.  The  following  modification 
to  Davidon* s  Variance  algorithm  provides  for  a  search  for  an  Improved 
estimate  for  x  at  each  iteration,  if  necessary,  using  several  computa¬ 
tions  of  f(x)  before  computing  g(x)  and  updating  the  estimate  of  G~^. 
Since  this  search  Introduces  a  major  change  in  the  algorithms,  the  form 
of  the  general  term  is  re-derived  as  follows. 

Assume  xn,  gjj,  Dq  *»  have  been  obtained  where  gj^  *  grad  fCXjj), 

DQ  is  positive  definite,  and  f(x)  =  tQ  a'x  +  1/2  x’Gx.  Pickanan  such 
that  fCxjj  -  OnDjjgjj)  <  f(Xjj)  using  *  1  whenever  possible,  otherwise 
Oji  is  reduced  (for  example  an  =  1/2,  1/4,...)  until  f  is  decreased. 

This  is  possible  since  [f  ’  (Xjj)  ,-Dngn]  =  -g’nDn8n<0'  Set  =  _ftr.DnV 
Xn+i  =  xn  +  °n  30 d  impute  gT1+1  =  grad  fO^)  =  gj,  +  Gon  and  yn  »  gn+1 
-  gH  =  Gon.  An  improved  estimate  to  G  of  the  form  Cn+^  =  Cn  +  unu^ 

is  sought  such  that  Cn+1an  ■=  Gon  =  y„  =  gn+1  -  gn-  As  before, 

rem  which  uR  =  gn+l  gn  °nan  f 


Cn+l°n  =  Cn°n  +  “n^n  =  Vi  "  gR  f 
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or 


*nfl  Z  (1  ~  °°)gn 

u’o 
u  n 


since  Cnon  =  -^(D^)  = 

For  clarity,  set  vn  *  g^,  -  (1  -  c^ig^  and  observe  Uq  has  the  for* 
°n 

2 

from  which  k  =  v'a_. 


Tn  .  Then  ~  *  *n  *  vn  =  yo 

k  k  Vn  Wk 


and  u'o  =  vngn  _ 

k  v’o„/k 


vn°n 


k. 


n”n* 


n  n 


Mow  ^nfl  *ay  be  =  C;J  +  vnvn  and  an  application  of  the 


v'o 

wnn 


rank-one  inversion  formula  yields  D^j.  3  ^n+1  *  D  “  PnvnvnDn 


v'o  +  v’D  v 
n  n  n  n  n 


Since  v;^  -  v^1[gnfl  -  gn  +  «ngQ] 

“  vnVn  +  ®nvn®n®n 
“  vnVn  "  Vn* 

D  v  v'D 

nay  be  simplified  to  *  Da  -  n  »  °  »  . 

vAVn 

This  form  of  the  algorithm  can  be  shown  to  have  the  same  properties 
as  Davidon's  form.  For  example  if  Dq  agrees  with  on  y,  i.e. ,. 

D^y  »  GT^y  *  a,  then  so  does  Dq+j.  This  is  demonstrated,  as  before, 
by  showing  v^Dny  =  0  as  follows: 

°nvn  "  V«n+1  “  *n  +  “nV  ’  V^n+l  “  «n>  *  ~n 

but  °n  “  *0+1  -  *n  "  G_1  <*Wl  "  V’  50  Vn  "  <Dn  "  ^Xgn+l  ’  V  * 
and  v^Dny  =  <gn+1  -  gn) '  (Dn  -  (T^y  =  0. 

This  form  of  the  algorithm  also  suffers  the  difficulty  of  Dn+^  not 
necessarily  being  positive  definite  at  each  iteration  for  an  arbitrary 
function.  This  difficulty  can  be  avoided  by  applying  Davidon's  method 


of  assuring  boundedness. 


un+l 


D  -  /8n+l^n8n+l 
n  v 


Dn+i  may  be  written  as 

) 


DnvnvnDn 


vnDnyn  ®n+l^n®n+l 
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or  D  .  =  D  + (A  -  l)DnvnvnDn. 

0+1  n  n  siuWi 

where  A_  =  1  -  8n+lPn8n+l  . 

v'D  v 
n  n7n 

Following  Davidon,  An  is  chosen  by  the  above  relation  whenever  possible, 
otherwise  such  that  a<.An<£  where  0<a<l<3<“  in  order  to  assure 
ax'  Dn+lx^  ’  Dnx^x’  Dnx 

for  any  x. 


3.  COMPLETE  RANK-ONE  ALGORITHM  FOR  FUNCTION  MINIMIZATION 

The  complete  algorithm  proceeds  as  follows: 

(1)  Initially  set  DQ  ■  I  and  pick  a,  f?  such  that  0<a<l<$. 

Choose  an  arbitrary  xD  and  compute  f(x0)  and  g(xQ)  ■  grad  f(xQ). 

(2)  Compute  f(xR  -  <V.DQ8n)  for  “n  =  1»  1^2»  '*•  until 

£<*n  “  ^nSn)  <  f<xn> 

(3)  Set  on  -  -anDngn  and  ^  -  x„  +  on 

Compute  gn+1  *  g(Xn+1)  -  grad  f(xn+1) 
and  set  yQ  -  gn+1  - 

vn  “  y„  +  Vn  “  *n+l  “  '  “n^n 

and  Yn  “  8n+lPn8n+l 
vn°nVn 

(4)  If  Yn>l  “  o  set  An  -  a, 
if  yn<l  -  g  set  An  =  g, 
otherwise  set  A  **  1  -  Yn* 

(5)  Set  Dn+1  =  Dn  4  <An  -  1)  DnvnvnDn 

®n+l^nSn+l 

(6)  Repeat  steps  (2)  through  (5)  until  selected  error  tolerances(s ) 
is  (are)  satisfied.  One  or  more  of  the  following  error  tests 
may  be  used: 
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I  I ®n+l 1 1  <  el*  8n+lDngn+l  <  V 
f(x.)  -  f(xn+1)  <  e3,  and  ||on|j  < 

It  is  interesting  to  note  that  any  algorithm  which  generates  re¬ 
cursive  estimates  to  the  inverse  of  the  Hessian  as  do  the  Fletcher- 
Powell  algorithm,  Davidon's  Variance  algorithm,  or  the  modified 
Davidon  algorithm  presented  above  can  be  expected  to  exhibit  numeri¬ 
cal  difficulties  whenever  the  Hessian  or  its  inverse  is  singular  at 
the  minimizing  point.  The  following  examples  illustrate  simple  problems 
with  this  property. 

Consider  the  problem  of  minimizing  f(x,y)  *  x2  y^.  The  first 
and  second  derivatives  are  Vf(x,y)  ■  / 2x\and 

W) 


V2f (x,y)  -  [2  0 

0  12y2  . 


Clearly  the  minimum  is  at  (x,y)  «*  (0,0)  where  Vf(x,y)  *  (0,0)'  but 
V2f(0, 


),0)  "J2  ®lis  singular.  Also, 

1°  °J 

^fCx.y)'1  - 


■p 

N5 

O 

_ 1 

1/2 

0* 

°  — , 

12y2J 

y-*0 

1 

0 

00 

Any  algorithm  which  attempts  to  estimate  V2f(x,y)“*  can  well  be  ex¬ 
pected  to  have  terms  which  tend  to  become  unbounded  as  the  solution  is 
approached. 

For  an  example  of  the  inverse  condition,  consider  minimizing 

f(x,y)  ■  x2  +  y^/3,  Then,  formally,  Vf(x,y)  =^2*^  1/3^ 

r  1 

V2f(x,y) 


and 


2  0 

0  A  "2/3 
9  y 


Clearly  the  minimum  is  at  (x,y)  =  (0,0)  and  in  this  case 
V2f(x,y)~1 


7TI 

fsj 

0 

_ 1 

1/2 

0 

0  |  y2/3 

—  m. 

y+0 

0 

0 
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In  this  case,  any  algoritha  which  attempts  to  generate  V^f(0,0)~^  can 
be  expected  to  become  singular. 

In  the  usual  formulation  of  the  Fletcher-Powell  algorithm  the 
matrix  Hr  is  tested  at  each  iteration  and  if  any  elements  become  too 
large  or  too  small  the  algorithm  is  re-started.  Thus,  in  the  worse 
case,  the  Fletcher-Powell  algorithm  would  de-generate  into  the  usual 
gradient  algorithm.  Davidon’s  algorithm,  both  the  original  and  as 
modified  here,  provide  assurance  (through  a, 3)  that  the  computations 
can  proceed  by  making  conservative  estimates  to  the  inverse  of  the 
Hessian. 

The  relative  speed  of  Davidon's  algorithm  compared  to  conjugate 
gradient  and  Fletcher-Powell  algorithms  was  reported  on  in  Reference 
11  .  Those  results  indicate  the  Davidon  method  is  the  superior 
algorithm  for  the  solution  of  most  of  the  test  problems  considered. 

In  the  few  cases  where  another  algorithm  was  found  superior  the 
differences  were  either  marginal  or,  as  in  a  few  cases,  the  Fletcher- 
Powell  algorithm  arrived  at  an  exact  solution. 
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SECTION  IV 

DIRECT  APPLICATION  OF  RANK-ONE 

Although  rank-one  methods  have  been  used  to  construct  algorithms 
for  minimizing  arbitrary  functions,  the  direct  application  of  rank-one 
to  minimizing  special  classes  of  functions  has  been  overlooked.  Con¬ 
sider,  for  example,  minimizing  the  function 

f(x)  =  1/ 2 1 J x 1 1  +  l/2k(a  -  m’x)2 

where  x,  ra£Rn  and  k,  a  are  scalars.  The  gradient  of  f  at  x  is 

Vf(x)  =  x  +  k(a  -  m’x)(-m). 


The  necessary  condition  for  f  to  be  minimum  at  x  =  x*  is  that  the 
gradient  (Vf(x*))  be  zero.  Furthermore,  since  m'x  is  a  scalar,  (m'x)m 
can  be  rewritten  as  m(m'x)  =  (tran’Jx,  where  mm’  is  the  outer  or  tensor 
product,  i.e.,  an  nxn  rank-one  matrix.  Use  this  fact  and  set  the 
gradient  to  zero  to  obtain 

x*  -  kam  +  kmm'x*  =  0 


(I  +  kmm')x*  =  kam  . 

/■ 

The  rank-one  matrix  inversion  technique  (Section  III)  is  applied  to 
obtain  directly 

x*  *  (I  +  kmm')  *kara 


-  (I 


kmm* 

1  +  km'm 


)kam 


ka 

B  «  *  i"  ■  i""  m 

1  t  km  m 


More  complex  forms  can  be  handled  with  little  difficulty.  Let  M 
be  a  m  x  n  matrix,  a  an  ra-vector  and  consider  minimizing 

f(x)  =  l/2||x||2n  +  l/2k| [ a  -  Mxj |2m 
=  l/2x'x  +  l/2k(a  Mx) ' (a  -  Mx). 
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The  gradient  of  this  function  is  given  by 

Vf  (x)  =  x  +  k(-M*)(a  -  Mx)  =  x  -  kM'a  +  kM'Hx. 
Again,  for  a  minimum  of  f  at  x* 

Vf(x*)  =  (I  +  kM’H)x*  -  kM'a  =  0 


or 


x*  =  (I  +  kM'M)-1kH'a 


th 


To  clarify  the  structure  of  M'M,  let  be  the  i—  row  of  M  represented 
as  a  column  vector,  m^  the  j.—  element  of  m^  as  well  as  the  i,j— 
element  of  M. 


Since  |  M’M  |  y  -  I  Vkj, 
bl 


and 


vt 


»  =  r™  n 


“kl 

"^2 
•  •  • 

mkn 


[*kl*  mk21 


\Al  "kl\2  *  ’  *  “kAli 

•  i 

mk2mkl  mk2\l  *  *  ’  VAn 


["Vkl  IVnk2  •**  raknmknJ 

it  follows  that 
m 

M'il  =  l  • 

k=l 

m 

Now  I  +  kM'M  may  be  represented  as  I  +  £  km  m'  and  its  inverse  may  be 

i=l  1 

computed  recursively  by  repeated  application  of  the  rank-one  inversion 
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■ethod.  The  question  of  the  existence  of  (I  +  leM'M)  *  is  answered  in 
Theorem  II  which  follows  shortly.  Formally,  the  method  of  computing 
(I  +  ltM*M)  *  is  derived  recursively  as  follows: 

Set  C  =  I  D  =  cf 1  =  I 

o  o  o 


S  “  1  + 


,-1 


D  - 
o 


kD  m.nJD 
olio 

1  +  kn'D  e. 
1  o  1 


Ci-1  +  tai"i 


-1 


'i-1 


1  +  kI*iDi-l,,i 


C  -  I  +  kM’M  D  -  C*1  -  (I  +  kM’M)-1 

n  mm 

Note  D  is  computed  in  precisely  m  steps  (m  *  row  order  of  M)  and  the  C*s 

B  1 

need  not  be  computed  since  can  be  computed  from  m^  and  k.  For 

this  problem,  x*  =  D  (fcM'a)  is  the  required  solution  satisfying  Vf(x*)  =  0. 

1.  APPLICATION  TO  AN  AIRCRAFT  WING-ROOT  BENDING  PROBLEM 

The  above  rank-one  method  was  used  to  determine  the  optimum  air¬ 
craft  wing  control  surface  deflections  required  to  minimize  a  specified 
penalty  function  {References  8,9).  The  objective  of  the  problem  was  to  reduce 
wing-root  bending  moments  through  active  control  of  trailing  edge  control 
surfaces  on  the  wing.  In  a  physical  application  such  reduced  bending 
moment  loads  could  lead  to  reduced  structural  requirements,  thus  reducing 
the  aircraft  weight  and  improving  aircraft  performance,  or  alternatively 
lead  to  an  expansion  of  the  aircraft  operational  envelope.  An  immediate 
consequence  of  the  control  surface  deflections  is  to  change  other  im¬ 
portant  aircraft  wing  characteristics,  principally  lift,  pitching  moment 
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and  drag,  these  also  effect  aircraft  perfomaece. 

For  nost  appl i cat ioas, increased  lift  at  a  gives  azgle-of-attack 
is  desirable  oecaase  as  increase  is  vertical  acceleration  can  be  obtained 
with  less  increase  in  wing  angle-of-attack.  This  leads  to  less  drag 
and  hence  sore  efficient  operation. 

Any  change  in  wing-generated  pitching  swat  would  require  changes 
in  aircraft  trim  to  generate  balancing  roccrents.  Although  trim  changes 
can  be  accomplished  automatically ,  this  adds  system  complexity  and 
interface  problems.  The  desire  to  avoid  complexity  is  motivated  by 
the  usually  valid  supposition  that  the  more  complex  system  is  less  re¬ 
liable. 

From  drag  considerations  and  because  linear  aerodynamic  theory 
was  used,  it  was  desired  to  restrict  the  nagnitude  of  control  surface 
deflections.  The  use  of  linear  aerodynamic  theory,  in  addition  to  being 
conceptually  and  computationally  simpler,  was  necessary  to  provide  the 
linear  system  description. 

The  mathematical  model  of  this  problem  was  formulated  as  follows: 
Changes  in  wing  lift,  pitching  moment  and  root  bending  can  be  repre¬ 
sented,  assuming  linear  aerodynamics,  as  a  linear  function  of  the  control 
surface  deflections;  therefore  let  6  =  (6^,  6^,  •••,  6^)  represent  the 
deflections  of  the  wing  control  surfaces,  =  m^6  the  change  in  wing 
lift,  &C„  *=  m'5  be  the  change  in  wing  pitching  moment,  £C  =  m'6  the 
change  in  wing-root  bending  moment.  The  following  cost  function  was 
formulated  whi-.h  when  minimized  would  tend  to  minimize  wing-root  bending 
moment  while  holding  changes  in  pitching  momeni  small,  maximize  the 
change  in  lift  and  maintain  reasonable  control  surface  deflections: 
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j(§)  -  +  V*S*>2  ’  S(4Cl)  +  1/2H*!S2 

where  the  weigbtlag  factors,  tj,  rrf  1^,  are  assented  bam  a  priori. 
The  solution  to  this  problem  is  readily  obtained  by  a  single  applica¬ 
tion  of  the  rank-one  method. 

Alt hoc gh  the  solutions  are  valid  the  results  are  difficult  to 
interpret  due  to  the  fact  that  each  solution  represents  a  different  wing- 
lift  condition.  The  problem  was  reformulated  at  a  constant  lift  condi¬ 
tion.  Specifically,  an  increment  in  angle-of-attack  was  incorporated 
to  allow  the  wing  to  generate  the  sane  lift  with  or  without  controls 
deflected.  The  aore  complete  representations  for  the  changes  in  wing 
characteristics  are: 

ACK  m  m2S  +  Sfaa’ 


ACEB  -  *36  +  C*Ba°' 


Thus  for  constant  lift  AC^  *  0  and  o  * - 

a 

Then  the  changes  in  pitching  -:~aent  and  wing-root  bending  at  constant 
lift  are: 

CHQ  _ 

AC^  =  m^6  -  - o' 6  ««  n'6 

Ct  1  2 


V/Rga 

ACbr  =  id's  -  -  m'5  =  m'6  . 

118  3  1  3 


Finally,  from  linear  aerodynamic  theory,  the  induced  drag  of  a 
wing  is  minimum  when  the  wing  pressure  distribution  is  elliptical. 
Therefore,  an  additional  terra  was  induced  to  represent  the  deviation  of 
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che  predicted  prcssarc  discri&acioe  from  elliptical-  Specifically,  tee 

spacsfise  distribmtioa  as  a  fssactioo  of  rt*  isenmlized  seacispaa  ts  was 

represented  as  C|((j)  *  C£  (tj)  +  c-  (cj)a  +  c-  (ss)<$  and  the  elliptical 

c  a  a 


-£  (tj)  =  —  /  1  *  CJZ. 

ellip  * 


distrifeetlca  as  c«  (tj)  =  — —  /  1  -  m*.  Following  tbe  above  procedure 
for  ccrrectisj  the  aogie-ef— attack  change. 


a  (n) 


ft  (n) 

ellip 


-  CgCn)  «  j^L  /  i  -  rs2  -  C£o  (rs)J  -  jc^Cd)  -  -q- - at’J 


Ibis  is  discretized  for  k  spaovise  locations  (ijj,  i=l,-  .  -,k)  then  a 
measure  of  the  spanvise  distribution  deviation  from  elliptical  is  given 
*>y  ll*o  “  Afill2*  »*ere 

2C 


IMi-  _i  A  -  ^  ^  faj) 

.  .  C1  ^ni) 

'*  'lj  '  %  - "13 


Lhx 


and  the  norn  is  in  R^. 


Finally,  the  following  cost  function  was  formulated  to  nininize 
wing-root  bending,  change  in  pitching  eoaent,  control  surface  deflections 
and  spam/ise  distribution  error: 

J(6)  =  K^AC^)  +  *2(ACm)2  +  1/2|  |6!  J2  +  X5|lX0  -  A«|||k 

JK. 

=  Kj¥26  +  K305^6)2  +  l/26'6  +  K_(A0-  A6)'(A0-A6). 

The  minimum  of  this  function  is  obtained  by  computing  the  gradient  and 
equating  the  gradient  to  zero.  The  gradient  is  given  by 
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?J(I)  *  +  4  +  t5H’)(i0-M) 

*  (I  +  e^3^3  +  Kjl’ -  KjA'l^ 
aed  the  rales-  of  5  for  vfsicEs  7J(<s)  *  0  is  given  fey 
«  =  (I  +  £^5^  +  -  g^). 

A 

The  jdcmriaj  $■  can  them  be  obtained  fey  fcs-1  iterations  of  tine  raak- 
ooe  method. 

In  both  of  the  above  formulations  of  the  problem  the  weighting 
factors,  Ij's  in  the  cost  function  are  presuned  known.  As  is  often  the 
case  in  this  type  problem,  this  is  generally  not  true.  However,  the 
solution  to  the  problem  for  a  given  set  of  K^'s  can  fee  computed  extremely 
rapidly  on  a  digital  computer  using  the  rank-one  nethod,  raking  a 
systematic  search  over  a  vide  variation  of  the  weighting  factors  a 
practical  approach  to  obtaining  3  realistic  solution.  For  example,  a 
program  written  for  the  CDC  6600  computer  to  solve  the  above  problen 
computes  the  basic  natrices  and  vectors  required  (A,»3’s,  etc.),  incre¬ 
ments  the  Kji's  in  a  systematic  fashion  and  computes  over  3,000  optitsal 
6's  in  approximately  10.0  seconds  central  processor,  CP,  tine.  With 
such  a  volume  of  data,  care  had  to  be  taken  in  the  search  pattern  and 
output  format  to  avoid  an  overwhelming  output  volume. 

2.  DERIVATION  OF  THE  DIRECT  RANK-ONE  METHOD 

In  the  preceding  development  only  the  solution  of  the  necessary 
condition  for  an  extremum,  V£(x*)  =  0,  was  considered.  For  the  simplest 
form,  a  linear  combination  of  quadratic  terms  with  positive  coefficients, 
if  an  extremum  exists  it  must  be  a  minimum.  The  addition  of  linear  terms 
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cskcs  ao  particular  difficulty  siace  the  ^natfrati:  terns  can  be  ex¬ 
pected  to  eventmally  dbattoate.  For  the  more  gtaeral  probleat  where  die 
weigfctiag  factors  are  arbitrary,  the  rack— coe  method  also  provides 
some  inform*  tice  cc  the  natore  of  the  solution. 


Fa>r  any  of  the  above  examples  the  gradient  is  represented  by  the 
general  font  ?f(x)  =  Cx  +  a,  were  C  is  a  matrix  and  a  a  vector.  The 
second  derivative  is  "2i  (x)  =  C.  The  following  Lemma  provides  informa¬ 
tion  on  the  positive  definiteness  of  C. 


II 


If  Cj  and  Da  *  C^1  are  positive  definite. 


U)  ■  a’D 

n  n  n  n 


1)11+1  Cn*l  *  (Cn  +  1  *  “a"  1  +  kn^Da»jj 


and  1  +  k*nDn*n>0»  then  C,^  aEU*  ®n+.i  are  positive  definite 


Proof 

cor  an  arbitrary  u,u'Dnu  >  0  therefore  the  ratio 


c'Via 

u'D  u 
n 


u 4  Don*  Du 
□  a  n  n 


u'D  u  1  +  ka'D  n 
“  n  n  n 


1  - 


ka'D  a 
non 


(u’D  a  )' 
n  n 


1+klaaDnton  <u  ’  Dnu>  (Qn  Vn> 


2 

By  the  Cauchy-Schvarz  inequality  (u'Dnnn)  <_  (u'D^u) (m^D^n^)  where 
equality  holds  if  and  only  if  u  is  a  scalar  multiple  of  m  .  Clea'ly 

« _  II 

u  Dn+lU 

the  ratio  _  =  1  if  u'Dnmn  =  0  and  this  ratio  is  furthest  from  one 

u'Dnu 

when  u  is  a  multiple  of  m  .  in  which  case 
r  n 
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*'v**iu  »  i  .  KPg? 

1  +  to'B« 

nan 

-  1 

1  +  kat’D  b  s 
ana 

Let  a  *  Min  J  1,  (1  + 

and  0  *  Max  { 1.(1  +  ko'D  st  )~*  ( 

non 


then  au*Dnu  <  n'D^n  <_  Su'FLjti 


"»°*  for  Cn+l  =  Cn  +  **n*h  * 


v*C 


»+l 


1  + 


kfrr’m^2 

v‘Cnv 


1 + kv»  ,(vy . — - 


Since  C  =  D  -1  and  C  and  D  are  positive  definite 
n  n  n  n 


-  (v'Cnv)KDnnn>' 


where  equality  holds  if  and  only  if  v  is 
the  ratio  v  ^n-<-lv  is  one  if  v’m  =  0  and 

-Tv 

multiple  of  in  which  case 


a  scalar  multiple  of  m^. 
furthest  from  one  when  v 


Thus 
is  a 


v ' Cn+lv 
v’Cnv 


1  +  kra'Dm  . 
n  n  n 


Let  y  =  Min  {  1,  1  +  kra'D  m  } 

non* 


1*0 
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amd  X  *  Mar  f  1,  1  +  I 

th«  Yv’Cv  <  v'C^v  £  X  v'C^, 

Thus  the  denominator,  1  +  km*D  ■  ,  provides  a  necessary  and  suf- 

nan 

ficient  condition  for  the  invertibiiity  of  Vi'  i.e. ,  the  existence  of 
by  Lemma  1  and  a  sufficient  condition  for  and  to  be 

positive  definite  by  Lemma  II. 

Before  proceeding  to  a  proof  of  the  rank-one  method .  the  following 
observations  are  made  to  simplify  the  notation.  For  the  problem  of 
minimizing 

£W  -  IMIJn  +  klHa  "  +  k2NXo  '  ^IIJp  » 

the  gradient  is  given  by 

Vf(x)  «2i+  2k1(-M*)(a  -  Mr) 

+  2k^{~A * ) (Xq-Ax) 

-  2(1  +  kjM'M  +  k^A’AJx  -  ZkjH’a  -  2k2*'*0. 

If  x*  minimizes  f  then  Vf(x*)  *  0  and  if  (I  +  kjM'M  +  k2A’A) 
exists,  then  x*  =  (I  +  kjM’M  +  k^' AJ^Ck^’a  +  k2A’Xo).  Furthermore, 
if  (I  +  hjM'M  +  k2A'A)  is  positive  definite  then  f  has  a  local  minimum 
at  x*. 

As  before  let  n^,  i  =  l,2,...,m,  be  the  i—  row  of  M  and  1^1=1, 
...,p,  be  the  i—  row  of  A.  Then  (I  +  k^M’M  +  k2A'A)  can  be  written  as 

H 

I  +  l 
i»l 

where  N  =  m  +  p 


hi 
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|ki  1“1*****“ 

Ikj  i  *  -,»fp 


i  -  1,...,. 


i  -=  *frl,...,«frp. 


Jefine  and  recursively  for  k  =*  0,1,..., N  as  follows: 


Co*1 


D  *=  c"1  -  I 
o  o 


°k  -  Cfc-1  +  T*kvkvk  dk  “  1  +  nkvkDk-lvk 


“  Vl  ~  Vl  vkvkDk-l, 


Iheoren  II 


With  C^.d^  and  as  defined  above 
H 

(a)  =  I  +  l  n^v* 

i-1 

(b)  so  long  as  d.  i*  0,  i  »  1,.. .  ,k 


D.^  ”  (T1  and  in  particular  for  k  «*  N, 


dn  -  S  -  <x  +  J'WP 

i“l 

(c>  if  >  0  for  all  k  =  1,2,...,N,  is 
positive  definite. 


Proof 


(a)  By  construction: 


C  o  I 
o  » 


C1  “  I  +  VjVJ, 
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e2  =  ci  +  n2v2v£ 


I  +  VlVl  +n2V2Vi 


1  +  l  ^ivivi . 

i-i 


a  -  i  +  [  niv1v^ 

K  1-1 


N 


S  "  I  +  IwpL 

i-i 


(b)  By  induction: 

Initially  CQ  «  I  at-  DQ  «  C0*  ■  1.  Suppose  the  assertion 
holds  for  k-1,  i.e.,  D^..^  **  and  d^  j*  0,  for  i»l,...,k-l.  Now  if 
dfc  J*  0,  exists  by  Lemma  I  and 

Hk 

‘VHc  “  (Ck-l  +  nkvkvk>  ^°k-l  _  —  ViWVi) 

‘  °k-l Vi"  ^  Ck-lI>k-lvkvkDk-l 

+  WkVl  "  —  vkvkDk-lvkvk°k-l . 

dk 

Since  *  I  and  is  a  scalar,  the  above  may  be  rewritten 

39  nk 

CkDk  "  1  “  d^  vkvk°k-l  +  nkvkvk°k-l 

n.2 
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■  1  "  vkvkDk-l  (1  ‘  dk  +  nkvk°k-lvk> 

■  It 

Where  a  1  +  is  used.  Thus  the  assertion  holds  for  k. 

Furthermore,  if  d^  f  0  for  k=l,...,N  then  =  Cfc*  for  k=l,...,N  and 
in  particular,  from  (a), 

N 

»  *1  Vi'l)'1 

i"l 

(c)  By  Induction: 

Clearly  DQ  “  C~*  ■  I  is  positive  definite. 

Suppose  the  assertion  holds  for  k-1,  i.e.,  is  positive  definite. 

If  dk  >  0  then  by  Lemma  II  the  ratio  u.  .PlcU...  lies  between  1  and 

u’Vlu 

d£*  >  0  for  arbitrary  u.  Hence,  u'D^u  is  positive  for  all  u  which  inr* 
is  positive  definite  and  the  assertion  holds  for  k.  If  dk  >  0  for 
all  k  -  1,...,N  then  Djj  is  positive  definite. 

It  should  be  noted  that  the  matrices  need  not  be  computed  since 
and  d^  can  be  computed  from  and  directly.  The  denominator 

d^  provides  a  convenient  check:  if  d^  is  zero  (or  sufficiently  small) 
then  is  singular  and  fails  to  exist  (D^  may  be  numerically  un- 
trac table) ;  if  d^  is  positive  at  every  iteration  then  the  resultant  ex¬ 
tremum  is  the  minimum  since  is  positive  definite.  However,  if  d^  is 
negative  for  one  or  more  iterations,  may  still  be  positive  definite 

since  Lemma  II  provides  only  a  sufficient  condition.  In  fact  it  may 
N 

happen  that  I  +  J  n  v  v'  is  invertible,  and  even  positive  definite, 
i“l  *  *  * 

however  an  intermediate  may  fail  t*'  exist  because  d^  is  zero.  In 


AF?DL-S?-7  2-77 


such  a  case,  a  reordering  of  the  steps  will  correct  the  difficulty  as 
shown  below. 

Suppose  all  steps  up  to  the  n^-»  n<N,  have  been  accomplished  and 

dQ  =  1  -  n^n^n-l^  =  0.  Then  Dq  fails  to  exist.  The  role  of  nnvnvp 

***  nn+lvn+lvn+l  ’>e  interchanged,  then  dn  =  1  +  nn+1v^+1Dn+1vn+1  is 

computed  and  presumed  non-zero.  Then 
—  ^n+l 

°n  “  Dn-1  "  ln  Dn-lvn+lv  .-1  Dn-1 
exists  and  is  the  inverse  of  Cr=  +  np+1vn+1v'+1# 


Now  the  w-rm  nnvnvp  is  again  considered 


d  *  1  +  0  v'D  v 
n+1  n  n  n  n 

nn+l 

*  1  +  *nvA<Dn-l  *  —  VlvA+lDn>vn 
dn 


t  T*nr|n+1  , 

1  +  n  v'D  ,v  -  - v  b  ,v  ,,v'  ,Dv 

n  n  n-1  n  d,^  n  11-1  n+1  n+1  n  n 


nnI1n+l  2 

-  ~y—  (vnDn-lvn+l>  > 
dn 


since  the  first  two  terms  form  dp  =  0.  In  general  dn+j  4  0  and  hence 

vrvi" 


dn+: 


DnvnvADn  exists* 


A  simple  numerical  example  is  presented  to  illustrate  this  special 
case.  Consider  the  problem  of  inverting  I  -  (q)(1,0)  +  (J)  (1,1). 
Proceeding  as  usual,  set  Dq  =  I  and  compute  d^  =  1  +  (-1)(1,0)(q)  = 
1-1  =  0. 


^5 
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Since  fails  to  exist,  the  role  of  (1,0)  and  (1,1)  are  inter¬ 
changed  and  d^  is  computed  as  d^  *  1  +  (1)(1,1)  (j)  ■  1  +  2  ■  3. 
Therefore,  =*  I  -  y(J)(l,l)  =  [*2/3  -1/3' 

-1/3  2/3 

J  * 

Now,  picking  up  the  first  term,  ^  ”  1  +  (-1) (1,0)D^(q)  ■  1  -  ^ 


(-1)  _  ! 

and  D2  *  D1  “  Yjr  Dl(0)(1*0)Di 


The  correctness  of  the  solution  is  readily  verified  as  the  inverse  of 


I  -  (q)(1,0)  +  (J) (1,1) 


b6 
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SECTION  W 
CONCLUSIONS 

The  main  topic  of  this  dissertation,  computational  methods  for  the 
solution  of  unconstrained  minimization  problems,  is  covered  in  three 
parts;  generalized  descent  algorithms,  a  rank-one  minimization  technique 
for  the  general  problem^and  a  direct  rank-one  method  for  a  special  class 
of  problems.  Although  each  topic  relates  to  unconstrained  minimization, 
the  scope  is  successively  decreased  as  the  specialization  is  increased. 

The  theorem  on  generalized  descent  algorithms.  Theorem  I,  demon¬ 
strates  the  essential  properties  of  a  descent  algorithm.  Here  a  descent 
algorithm  is  defined  to  be  a  computational  method  in  which,  at  each  iter¬ 
ation,  a  descent  direction  is  generated  and  a  single-dimensional  search 
is  conducted  for  a  minimum.  Since  the  setting  is  highly  abstract,  addi¬ 
tional  specialization  is  included  and  the  theorem  is  applied  to  the  three 
descent  algorithms  in  common  ,use:  the  gradient  method,  conjugate  gradients 
and  the  Fletcher-Pdwell  method.  The  essential  property  of  descent 

I 

algorithms,  choosing  a  ’escent  direction  and  the  search  for  a  minimum  at 
each  iteration,  is  sufficient  to  cause  the  sequence  of  function  values 
to  monotonically  decrease.  An  additional  property  is  required  to  assure 
the  derivative  of  t  e  function  goes  to  zero. 

The  descent  direction  generated  by  the  algorithm  must  be  strict  in 
the  sense  that  convergence  of  the  inner  product  [f '  (x) , <f> (x)  ]  .  to  zero 
must  imply  convergence  of  the  derivative  to  zero.  This  condition  is  made 
rigorous  in  condition  (ii)  required  of  the  function  This  condition 
and  the  uniform  continuity  of  the  derivative  f'  are  used  to  demonstrate 
the  derivative  must  converge  to  zero. 


^7 
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Although  Theorem  1  demonstrates  the  common  properties  of  descent 
algorithms ,  It  fails  to  provide  any  information  on  one  very  important 
property  of  ail  such  algorithms,  convergence  rate.  Convergence  rate 
involves  not  only  the  amount  of  improvement  at  each  iteration  but  also, 
for  very  practical  reasons,  the  computational  time  required  to  accomplish 
the  iterations. 

One  of  the  common  properties  of  descent  algorithms,  the  single¬ 
dimensional  search  can  also  be  a  drawback  in  that  this  is  generally  the 
most  time  consuming  step  of  each  iteration.  Therefore  an  algorithm 
which  minimizes  a  function  without  requiring  repeated  single-dimensional 
searches  for  a  minimum  might  be  superior  to  any  descent  method.  One 
such  method,  Davidon's  Variance  Algorithm  (Reference  13),  has  been  shown  to  be,  in 
many  cases,  superior  to  the  three  common  descent  algorithms  (Reference  11). 
Davidon's  method  is  discussed  in  Section  III. 

A  new  derivation  of  Davidon's  method  is  presented.  This  deriva¬ 
tion  provides  a  clearer  insight  into  the  structure  of  the  algorithm  by 
considering  the  algorithm  as  repeated  applications  of  the  rank-one  matrix 
inversion  technique.  Although  the  rank-one  method  would  be  exact  when 
applied  to  a  quadratic  fora  and  would  require  no  single-dimensional 
search,  certain  added  computational  precautions  must  be  included  if  the 
method  is  to  be  used  to  minimize  an  arbitrary  function.  One  of  these  is 
provided  by  Davidon  and  assures  the  matrix  approximation  to  the  inverse 
of  the  Hessian  is  well  behaved.  Another  precaution  provided  by  this 
author  assures  only  "good"  estimates  of  t.ie  solution  are  used  to  update 
the  approximation  of  the  inverse  of  the  Hessian.  Although  a  linear 
search  is  involved,  only  the  function  value  is  computed  and  the  function 
need  only  be  decreased  not  minimized.  As  a  result,  for  each  iteration 
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both  the  number  of  test  points  is  reduced  and  the  complexity  of  compu¬ 
tations  at  each  test  point  is  reduced. 

An  extension  of  the  Davidon  algorithm  to  the  problem  of  minimiz¬ 
ing  an  arbitrary  function  defined  on  an  infinite  dimensioned  space  such 
as  a  functior  space  is  highly  interesting.  In  the  finite  dimensional 
case  the  outer  product  of  vectors,  uv* ,  can  be  represented  as  a  natrix. 

If  u  and  v  belong  to  a  function  space,  the  representation  of  uv*  is  not 
entirely  obvious.  However,  since  most  problems  require  digital  computa¬ 
tion,  and  the  representation  of  a  function  is  necessarily  discretized, 
the  function  space  may  be  considered  as  Is  with  H  very  large. 

The  third  area  considered  in  this  dissertation  applies  the  rank- 
one  method  to  a  special  class  of  problems.  Although  the  class  of  prob¬ 
lems  to  which  the  method  applies  is  specialized,  it  is  not  uncommon. 

In  fact  many  preliminary  engineering  problems  are  in  this  class  where 
the  cost  function  is  a  combination  of  linear  and  quadratic  terms  each 
weighted  by  a  penalty  factor  which  is  constant  but  unknown.  The  rank- 
one  method  provides  a  rapid  solution  to  the  problem.  The  method  also 
provides  necessary  and  sufficient  tests  for  the  existence  of  an  extre¬ 
mum  and  a  sufficient  test  for  the  solution  to  be  a  minimum. 

Extension  of  the  direct  rank-one  method  to  problems  in  a  function 
space  is  of  continuing  interest  to  the  author.  As  noted  for  the  ex¬ 
tension  of  Davidon' s  method,  care  must  be  observed  in  the  interpreta¬ 
tion  of  the  outer  product.  Since  the  number  of  iterations  of  rank-one 
required  to  obtain  the  solution  is  determined  by  the  structure  of  the 
cost  function  to  be  minimized,  not  the  dimension  of  the  underlying  vector 
space,  application  of  the  method  in  a  function  sp.ace  has  possibilities. 
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To  date  the  class  of  problems  in  a  function  space  which  nap  be  solved 
by  the  rank-one  method  is  very  restricted  and  the  set  of  known  appli¬ 
cations  is  empty.  Finally,  the  problen  of  representing  functions  on  a 
digital  computer  still  exists. 

To  sumsarize,  the  following  extensions  appear  to  be  promising 
areas  for  further  investigation.  Application  of  Theorem  1  :o  other 
descent  algorithms-  Application  of  Dari don's  rank-one  aethod  to  a  large 
sanple  of  test  problems  to  determine  in  a  practical  application  its 
computational  speed.  Application  of  Davi don's  aethod  in  a  function  space 
night  be  fruitful  if  first  the  direct  rank-one  aethod  can  be  successfully 
applied. 


50 


i3r?m^z-j2~n 


agTOPwitras 


I-  F.  S-  3eckmac,  "Sfce  Salntice.  ©f  linear  EqzsZiczs  by  -.be  Conjugate 
Gradient  Method*,  Mshhiemaftical  Methods  for  Digital  Oan^cters, 

A-  Salstoit  and  5.  S.  Vilf  (Sds.),  John  tfiley  4  Sons,  Inc-,  few 
Tori,  S-  T-,  19&>- 

2.  if.  C  Barite:,  Variable  Metric  Methcd  for  Kiz.im.izs.ticc.,  Argocme 
S2.ti0tr.yJ  Lab- ,  Seoort  599©  (3er.  J,  1®59- 

3-  V.  C.  Barite:,  “Variance  Algorithm.  for  JtirEiorizaticc",  Computer  J-, 
10  a 965),  to.  lOS-llO. 

fc.  H.  Fletcher  sod  X.  J.  B.  Powell,  "A  Eapidiy  Chare rgent  Descent 
Method  for  Minimization”,  Ccarruter  J. ,  6  (l$63),  pp-  io3-!-S. 

5-  B.  Fletcher  and  C.  X.  Seeres,  “Function  Minimization  by  Conjugate 
Gradients*,  Ccanrafcer  J. ,  7  (19&)*  PP-  1S9-15'- 

6.  A-  A.  Goldstein,  "On  Steepest  Descent”,  J.  SIAX  Control,  Ser.  A, 
Vol.  3,  5c-  1,  1965,  pp-  1*7-1;:. 

7-  A-  A.  Goldstein,  “Minimizing  Functionals  or.  Sorted  Linear  Spaces". 
J.  SIAX  Control,  Yol.  U ,  So.  1,  2966,  pc.  81-89- 

8.  J.  2.  Jenkins,  D.  C.  Sckholat ,  ;md  3.  T.  Kujavski,  "An  Assess¬ 
ment  of  the  Interfacing  Problem;  with  CCY  Design  Concepts”,  AIAA 
2nd  Aircraft  Design  and  Opera .20ns  Meeting,  July  1970,  AIAA  Paper 
Ho.  70-926. 

9.  B.  T.  Khjawski,  J.  2.  Jenkins,  and  D.  C.  Eckholdt,  "Longitudinal 
Analysis  of  Two  CCV  Design  Concepts",  AIAA  3rd  Aircraft  Design  and 
Operations  Meeting,  July  1971,  AIAA  Paper  Ho.  71-786. 

10.  L.  S.  Lasdon,  S.  II.  Hitter,  and  A.  D.  Warren,  "I5ie  Conjugate 
Gradient  Method  for  Optimal  Control  Problems",  IEEE  Trans,  on  Auto¬ 
matic  Coitrol,  Vol.  AC-12,  Nr.  2,  April  1967,  PP-  132-138. 

11.  Craig  E.  Miller,  A  Computational  Comparison  of  Gradient  Minimiza¬ 
tion  Algorithms,  Masters  Thesis.  Wright-Patterson  Air  Force  Base, 
Ohio:  Air  Force  Institute  of  Technology,  March  1971- 

12.  M.  J.  D.  Powell,  "A  Survey  of  Numerical  Methods  for  Unconstrained 
Optimization",  SIAM  Review,  Vol.  12,  Nr.  1,  (1970)  pp.  79-97- 

13-  A.  E.  Taylor,  Advanced  Calculus,  Ginn  and  Company,  Boston,  Mass., 
1955. 

lU.  M.  M.  Vainberg,  Variational  Methods  for  the  Study  of  Nonlinear 
Operators ,  Holden-Day,  Inc.,  196*. 

15.  H.  S.  Wilf,  "Matrix  Inversion  by  the  Annihilation  of  Rank", 

J.  SIAM,  Vol.  7,  Nr.  2,  June  1959,  PP-  1*9-151- 


51 


